import sys
import cplex
import numpy as np
from cplex.exceptions import CplexError

ca = cplex.Cplex()
ca.read(sys.argv[1])

#copy ca to c for self processing
c = cplex.Cplex(ca)

c.parameters.output.writelevel.set(0)

# record index of binary variables
int_var_idx = []
con_var_idx = []

for i in range(c.variables.get_num()):
    if c.variables.get_types(i) != 'C':
        int_var_idx.append(i)
    else:
        con_var_idx.append(i)
        
# change the type from MILP into LP (0:LP, 1:MILP)
c.set_problem_type(type=0)

ub = c.variables.get_upper_bounds()
lb = c.variables.get_lower_bounds()

sense = c.objective.get_sense()

# three vectors to track the progress of tree search process
# use three list to store the information

P = []   # P: path vector
S = []   # S: sign vector (the sign indicate the bound value)
V = []   # V: value vector

s = 0
v    = 0.0
v_up = 0.0
v_lo = 0.0
v_fr = 0.0

branch_var = 0
branch_dir = 0

# consider min and max problem
# declare global z and node z
if sense == -1:  #maximization
    z_global = -10000000
else: #minimization
    z_global = 10000000

z_node = 0

# c.variables.set_lower_bounds to fix values at 0 or 1
# c.solution.get_status() to know the status for logical control  
# 1. optimal
# 2. unbounded
# 3. infeasible 

# convert output from cplex into numpy ndarray format
#       x = np.array(x)      
#
# Branching variable selection:
# For branching on the most fractional variable
# Test whether solution of LP is integer feasible
# if yes, fathom the node by integrality
# if false, we should trace the index of most fractional
# variable
#
#      s = argmax(np.minimum( x - np.floor(x) , np.ceil(x) - x) )
#
# then branch up by modifying upper bound of x_s
#
#      c.variables.set_upper_bounds(s, 1)
#
# or branch down by modifying lower bound of x_s
#
#      c.variables.set_lower_bounds(s, 0)

def backtrack(c, P, S, V, lb, ub):
    print 'Starting backtrack'
    print 'P', P, 'S', S, 'V', V
    idx = 0

    if len(P) == 0:
        print 'Already at the root'
        return []

    # check the sign last elements of P
    while P[-1] > 0:
        # skip through all fathomed nodes
        idx = abs(P[-1]) - 1
        
        if S[-1] == -1:
            c.variables.set_upper_bounds( idx, ub[idx] ) #relax the upper bound
        else:
            c.variables.set_lower_bounds( idx, lb[idx] ) #relax the lower bound
        P.pop()
        S.pop()
        V.pop()
        print 'backtrack one level'
        print 'P', P, 'S', S, 'V', V
        if len(P) == 0:
            break

    # if back to the root, return 
    if len(P) == 0:
        return []

    # now the idx is changed 
    idx = abs(P[-1]) - 1
    # then change the sign of S
    if S[-1] == -1:
        c.variables.set_upper_bounds( idx, ub[idx] )
        c.variables.set_lower_bounds( idx, V[-1] + 1 )
        V[-1] = V[-1] + 1
    else:
        c.variables.set_lower_bounds( idx, lb[idx] )
        c.variables.set_upper_bounds( idx, V[-1] - 1 )
        V[-1] = V[-1] - 1

    P[-1] = -P[-1]
    S[-1] = -S[-1]

    print 'Finish backtracking'
    print 'P', P, 'S', S, 'V', V
    
    return []
#End of backtrack function


iter = 0

c.parameters.simplex.display.set(0)

while True:
    # be ready to solve the node LP
    iter = iter + 1

    c.solve()
    print 
    print 'Iteratation', iter
    print 'status', c.solution.get_status_string(), c.solution.get_status()
    #print 'current upper bound', c.variables.get_lower_bounds()
    #print 'cuurent lower bound', c.variables.get_upper_bounds()

    # check the solution status
    if c.solution.get_status() == 3:
        # fathom by infeasibility
        # update P and S and add modify c
        print 'fathom by infeasiblity'
        backtrack(c, P, S, V, lb, ub)
    else:
        # get solution at node LP
        x = np.array( c.solution.get_values() )
        
        if not (x - np.floor(x))[int_var_idx].any():
            # integer feasible then update the node optimal value
            print 'fathom by optimality'
            z_node = c.solution.get_objective_value()
            if sense == -1: #maximization
                z_global = max(z_node, z_global)
            else: #minimization
                z_global = min(z_node, z_global)

            # fathom by integrality and then backtrack
            backtrack(c, P, S, V, lb, ub)

        else:
            # integer infeasible then check the bounds
            z_node = c.solution.get_objective_value()

            if z_node <= z_global and sense == -1:  #maximization
                # fathom by bounds and backtrack
                print 'fathom by bounds'
                backtrack(c, P, S, V, lb, ub)
            elif z_node >= z_global and sense == 1: #minimization    
                print 'fathom by bounds'
                backtrack(c, P, S, V, lb, ub)

            else:                
                # select branch variable
                
                x_int_frac_min = np.minimum( x - np.floor(x) , np.ceil(x) - x)
                x_int_frac_min[con_var_idx] = 0

                s = np.argmax( x_int_frac_min )

                # v = c.solution.get_values(s)
                v = x[s]
                v_up = np.ceil(v)
                v_lo = np.floor(v)
                v_fr = v - v_lo
            
                ### update P and S and modify c and go the next node always go left
                P.append(-s - 1)

                # print 's', s
                branch_var = s
                if v_fr <=  0.5:
                    # go left
                    S.append(-1)
                    V.append(v_lo)
                    c.variables.set_upper_bounds(s, v_lo)
                    branch_dir = -1
                else:
                    #go right
                    S.append(1)
                    V.append(v_up)
                    c.variables.set_lower_bounds(s, v_up)
                    branch_dir = +1

                    # print c.variables.get_lower_bounds()

                print 'branching variable:', branch_var, 'branching direction', branch_dir
    print '----------------'
    print 'z_node'  , z_node
    print 'z_global', z_global

    if len(P) == 0:
        break
#End of while loop

print 
print 'Final solution:'
print 'z_node', z_node
print 'z_global', z_global


# post processing    
#sol = c.solution

#print   
# print out problem info

# solution.get_status() returns an integer code
#print "Solution status = " , sol.get_status(), ":",
#print sol.status[sol.get_status()]
#print "Solutions: ", sol.get_values()
#print "Objective values: " , sol.get_objective_value()

